Degradation statistics

baRulho:quantifying habitat-induced degradation of (animal) acoustic signals

Author

Marcelo Araya-Salas

Published

July 3, 2023

Source code, data and annotation protocol found at https://github.com/maRce10/barulho_paper

Purposes

  • Estimate degradation metrics on re-recorded signals from playback experiment at Bosque de Tlalpan, Mexico City, 2019

  • Run statistical analyses

 

Load packages

Code
source("https://raw.githubusercontent.com/maRce10/sketchy/main/R/load_packages.R")

# install/ load packages
load_packages(packages = c("remotes", "kableExtra", "knitr", "formatR",
    "rprojroot", "maRce10/warbleR", "ggplot2", "tidyr", "viridis",
    "corrplot", "brms", "ggdist", "cowplot", "cmdstanr", "maRce10/brmsish",
    "emmeans", "ggsignif"))

my.viridis <- function(...) viridis(alpha = 0.5, begin = 0.3, end = 0.7,
    ...)

source("~/Dropbox/R_package_testing/brmsish/R/extended_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")
Code
degrad_df <- read.csv("./data/processed/tlalpan_degradation_metrics_v01.csv")

degrad_params <- c("blur.ratio", "spectrum.blur.ratio", "envelope.correlation",
    "excess.attenuation", "signal.to.noise.ratio", "cross.correlation",
    "tail.to.signal.ratio", "tail.to.noise.ratio", "spectrum.correlation")

comp.cases <- complete.cases(degrad_df[, names(degrad_df) %in% degrad_params])

pca <- prcomp(degrad_df[comp.cases, names(degrad_df) %in% degrad_params],
    scale. = TRUE)

# add to data
degrad_df$pc1 <- NA
degrad_df$pc1[comp.cases] <- pca$x[, 1]
degrad_df$pc1.1m.rate <- degrad_df$distance
# plot rotation values by PC
pca_rot <- as.data.frame(pca$rotation[, 1:4])

pca_rot_stck <- stack(pca_rot)

pca_rot_stck$variable <- rownames(pca_rot)
pca_rot_stck$Sign <- ifelse(pca_rot_stck$values > 0, "Positive", "Negative")
pca_rot_stck$rotation <- abs(pca_rot_stck$values)

ggplot(pca_rot_stck, aes(x = variable, y = rotation, fill = Sign)) +
    geom_col() + coord_flip() + scale_fill_viridis_d(alpha = 0.7,
    begin = 0.3, end = 0.8) + facet_wrap(~ind) + theme_classic()

1 Correlation among metrics

Raw metrics:

Code
cormat <- cor(degrad_df[, degrad_params], use = "pairwise.complete.obs")

rownames(cormat) <- colnames(cormat) <- degrad_params

cols_corr <- colorRampPalette(c("white", "white", viridis(4, direction = -1)))(10)

cp <- corrplot.mixed(cormat, tl.cex = 0.7, upper.col = cols_corr,
    lower.col = cols_corr, order = "hclust", lower = "number", upper = "ellipse",
    tl.col = "black")

Code
# sort parameters as in clusters for cross correlation
degrad_params <- degrad_params[match(rownames(cp$corr), degrad_params)]

2 Data description

  • 20 frequencies
  • 3 locations
  • 11520 test sounds
  • 160 sound treatment combinations
  • Sample sizes per location, transect and signal type
Code
agg <- aggregate(cbind(sound.id, treatment.replicates) ~ location +
    habitat.structure + distance, degrad_df, function(x) length(unique(x)))

agg$replicates <- agg$sound.id/agg$treatment.replicates

df1 <- knitr::kable(agg, row.names = FALSE, escape = FALSE, format = "html")

kable_styling(df1, bootstrap_options = c("striped", "hover", "condensed",
    "responsive"), full_width = FALSE, font_size = 15)
location habitat.structure distance sound.id treatment.replicates replicates
1 closed 10 480 160 3
2 closed 10 480 160 3
3 closed 10 480 160 3
1 open 10 480 160 3
2 open 10 480 160 3
3 open 10 480 160 3
1 closed 30 480 160 3
2 closed 30 480 160 3
3 closed 30 480 160 3
1 open 30 480 160 3
2 open 30 480 160 3
3 open 30 480 160 3
1 closed 65 480 160 3
2 closed 65 480 160 3
3 closed 65 480 160 3
1 open 65 480 160 3
2 open 65 480 160 3
3 open 65 480 160 3
1 closed 100 480 160 3
2 closed 100 480 160 3
3 closed 100 480 160 3
1 open 100 480 160 3
2 open 100 480 160 3
3 open 100 480 160 3

3 Statistical analysis

Code
iter <- 10000
chains <- 4
priors <- c(prior(normal(0, 6), class = "b"), prior(cauchy(0, 4),
    class = "sd"))


# set frequency to mean-centered
degrad_df$frequency <- degrad_df$frequency - mean(degrad_df$frequency)


# set base level for factors
degrad_df$habitat.structure <- factor(degrad_df$habitat.structure,
    levels = c("open", "closed"))
degrad_df$frequency.modulation <- factor(degrad_df$frequency.modulation,
    levels = c("no_fm", "fm"))
degrad_df$amplitude.modulation <- factor(degrad_df$amplitude.modulation,
    levels = c("no_am", "am"))
degrad_df$duration <- factor(degrad_df$duration, levels = c("short",
    "long"))
degrad_df$location <- as.factor(degrad_df$location)

degrad_df$distance_f <- paste0(degrad_df$distance, "m")
degrad_df$distance_f <- factor(degrad_df$distance_f, levels = c("10m",
    "30m", "65m", "100m"), ordered = TRUE)


set.seed(123)

cmdstanr::set_cmdstan_path("~/Documentos/cmdstan/")

# to run within-chain parallelization
mod_pc1 <- brm(formula = pc1 ~ frequency * habitat.structure + frequency.modulation *
    habitat.structure + amplitude.modulation * habitat.structure +
    duration * habitat.structure + mo(distance_f) + (1 | location) +
    (1 | treatment.replicates), data = degrad_df, prior = priors,
    iter = iter, chains = chains, cores = chains, control = list(adapt_delta = 0.99,
        max_treedepth = 15), file = "./data/processed/regression_model_pc1.RDS",
    file_refit = "always", backend = "cmdstanr", threads = threading(10))

mod_blurratio <- brm(formula = blur.ratio ~ frequency * habitat.structure +
    frequency.modulation * habitat.structure + amplitude.modulation *
    habitat.structure + duration * habitat.structure + mo(distance_f) +
    (1 | location) + (1 | treatment.replicates), data = degrad_df,
    prior = priors, iter = iter, chains = chains, cores = chains,
    control = list(adapt_delta = 0.99, max_treedepth = 15), file = "./data/processed/regression_model_blur_ratio.RDS",
    file_refit = "always", backend = "cmdstanr", threads = threading(10))

3.1 PC1 degradation

Code
effects <- c("habitat structure", "frequency modulation", "amplitude modulation",
    "frequency", "duration", "frequency:habitat structure", "habitat structure:duration",
    "habitat structure:amplitude modulation", "habitat structure:frequency modulation")

mod <- readRDS("./data/processed/regression_model_pc1.RDS")

ggs <- extended_summary(mod, n.posterior = 1000, fill = viridis(10)[7],
    trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE,
    gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam",
        "frequency.modulationfm", "durationlong"), gsub.replacement = c("",
        "habitat structure", "amplitude modulation", "frequency modulation",
        "duration"), effects = effects, print.name = FALSE, trace = TRUE,
    return = TRUE)

3.1.1 Conditional plots

Code
bs <- 10
ts <- 4

ggf <- plot(conditional_effects(mod, "frequency:habitat.structure"),
    plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Open",
    "Closed")) + theme_classic(base_size = bs) + labs(x = "Frequency",
    y = "Degradation (PC1)", color = "Habitat") + scale_fill_discrete(guide = "none") +
    theme(legend.position = c(0.2, 0.95)) + scale_x_continuous(breaks = seq(-5,
    5, 2.5), labels = round(seq(-5, 5, 2.5) + mean(read.csv("./data/processed/tlalpan_degradation_metrics_v01.csv")$frequency),
    1)) + geom_text(x = 0, y = -0.5, label = "*", size = ts)

ggfm <- plot(conditional_effects(mod, "habitat.structure:frequency.modulation"),
    plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Tonal",
    "FM")) + theme_classic(base_size = bs) + labs(x = "Habitat", y = "",
    color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
    0.5)) + ylim(c(-5, 1.5))

diff <- 0.1

# add contrasts
for (i in 1:2) ggfm <- ggfm + geom_signif(y_position = c(0), xmin = c(i -
    diff), xmax = c(i + diff), annotation = c("*"), tip_length = 0.02,
    textsize = ts, size = 0.4, color = "gray40")

ggfm <- ggfm + geom_signif(y_position = c(1), xmin = 1, xmax = 2,
    annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4,
    color = "gray40")

ggam <- plot(conditional_effects(mod, "habitat.structure:amplitude.modulation"),
    plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Flat",
    "AM")) + theme_classic(base_size = bs) + labs(x = "Habitat", y = "Degradation (PC1)",
    color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
    0.5)) + ylim(c(-5, 1.5))

# add contrasts
for (i in 1:2) ggam <- ggam + geom_signif(y_position = c(0), xmin = c(i -
    diff), xmax = c(i + diff), annotation = c("*"), tip_length = 0.02,
    textsize = ts, size = 0.4, color = "gray40")

ggam <- ggam + geom_signif(y_position = c(1), xmin = 1, xmax = 2,
    annotation = c("ns"), tip_length = 0.02, textsize = ts, size = 0.4,
    color = "gray40")


ggdur <- plot(conditional_effects(mod, "habitat.structure:duration"),
    plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Short",
    "Long")) + theme_classic(base_size = bs) + labs(x = "Habitat",
    y = "", color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
    0.5)) + ylim(c(-5, 1.5))

# add contrasts
for (i in 1:2) ggdur <- ggdur + geom_signif(y_position = c(0), xmin = c(i -
    diff), xmax = c(i + diff), annotation = c("ns"), tip_length = 0.02,
    textsize = ts, size = 0.4, color = "gray40")

ggdur <- ggdur + geom_signif(y_position = c(1), xmin = 1, xmax = 2,
    annotation = c("ns"), tip_length = 0.02, textsize = 4, size = 0.4,
    color = "gray40")


plot_grid(ggf, ggfm, ggam, ggdur)

3.2 Blur ratio

Code
mod <- readRDS("./data/processed/regression_model_blur_ratio.RDS")

extended_summary(fit = mod, n.posterior = 1000, fill = viridis(10)[7],
    trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE,
    gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam",
        "frequency.modulationfm", "durationlong"), gsub.replacement = c("",
        "habitat structure", "amplitude modulation", "frequency modulation",
        "duration"), effects = effects, print.name = FALSE)
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 6) Intercept-student_t(3, 0.1, 2.5) sd-cauchy(0, 4) sigma-student_t(3, 0, 2.5) simo-dirichlet(1) blur.ratio ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates) 10000 4 1 5000 187 0 1772.987 2746.438 497531395
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
habitat structure 0.029 0.026 0.033 1.001 14071.545 14362.080
frequency modulation 0.065 0.058 0.072 1.001 2201.419 3919.558
amplitude modulation 0.024 0.017 0.031 1.002 2032.491 4248.172
frequency 0.001 0.000 0.002 1.002 2355.250 4752.024
duration 0.001 -0.006 0.008 1.002 1772.987 3901.099
frequency:habitat structure 0.003 0.002 0.003 1.001 7870.588 2746.438
habitat structure:duration 0.005 0.001 0.009 1 21697.716 14527.390
habitat structure:amplitude modulation 0.002 -0.001 0.006 1.001 19309.028 14793.951
habitat structure:frequency modulation -0.020 -0.024 -0.017 1 22550.716 13648.971

3.2.1 Conditional plots

Code
bs <- 10
ts <- 4

ggf <- plot(conditional_effects(mod, "frequency:habitat.structure"),
    plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Open",
    "Closed")) + theme_classic(base_size = bs) + labs(x = "Frequency",
    y = "Degradation (PC1)", color = "Habitat") + scale_fill_discrete(guide = "none") +
    theme(legend.position = c(0.2, 0.95)) + scale_x_continuous(breaks = seq(-5,
    5, 2.5), labels = round(seq(-5, 5, 2.5) + mean(read.csv("./data/processed/tlalpan_degradation_metrics_v01.csv")$frequency),
    1)) + geom_text(x = 0, y = 0.16, label = "ns", size = ts)

ylim <- c(-0.1, 0.5)

ggfm <- plot(conditional_effects(mod, "habitat.structure:frequency.modulation"),
    plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Tonal",
    "FM")) + theme_classic(base_size = bs) + labs(x = "Habitat", y = "",
    color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
    0.5)) + ylim(ylim)

diff <- 0.1


# add contrasts
for (i in 1:2) ggfm <- ggfm + geom_signif(y_position = c(0.25), xmin = c(i -
    diff), xmax = c(i + diff), annotation = c("*"), tip_length = 0.02,
    textsize = ts, size = 0.4, color = "gray40")

ggfm <- ggfm + geom_signif(y_position = c(0.3), xmin = 1, xmax = 2,
    annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4,
    color = "gray40")

ggam <- plot(conditional_effects(mod, "habitat.structure:amplitude.modulation"),
    plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Flat",
    "AM")) + theme_classic(base_size = bs) + labs(x = "Habitat", y = "Degradation (PC1)",
    color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
    0.5)) + ylim(ylim)

# add contrasts
for (i in 1:2) ggam <- ggam + geom_signif(y_position = c(0.2), xmin = c(i -
    diff), xmax = c(i + diff), annotation = c("*"), tip_length = 0.02,
    textsize = ts, size = 0.4, color = "gray40")

ggam <- ggam + geom_signif(y_position = c(0.3), xmin = 1, xmax = 2,
    annotation = c("ns"), tip_length = 0.02, textsize = ts, size = 0.4,
    color = "gray40")


ggdur <- plot(conditional_effects(mod, "habitat.structure:duration"),
    plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Short",
    "Long")) + theme_classic(base_size = bs) + labs(x = "Habitat",
    y = "", color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
    0.5)) + ylim(ylim)

# add contrasts
for (i in 1:2) ggdur <- ggdur + geom_signif(y_position = c(0.2), xmin = c(i -
    diff), xmax = c(i + diff), annotation = c("ns"), tip_length = 0.02,
    textsize = ts, size = 0.4, color = "gray40")

ggdur <- ggdur + geom_signif(y_position = c(0.3), xmin = 1, xmax = 2,
    annotation = c("ns"), tip_length = 0.02, textsize = 4, size = 0.4,
    color = "gray40")


plot_grid(ggf, ggfm, ggam, ggdur)

3.3 Excess attenuation

Code
mod <- readRDS("./data/processed/regression_model_excess_attenuation.RDS")

extended_summary(fit = mod, n.posterior = 1000, fill = viridis(10)[7],
    trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE,
    gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam",
        "frequency.modulationfm", "durationlong"), gsub.replacement = c("",
        "habitat structure", "amplitude modulation", "frequency modulation",
        "duration"), effects = effects, print.name = FALSE)
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 6) Intercept-student_t(3, 83.8, 34) sd-cauchy(0, 4) sigma-student_t(3, 0, 34) simo-dirichlet(1) excess.attenuation ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates) 10000 4 1 5000 9 0 2625.42 5199.726 865228054
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
habitat structure 9.891 9.052 10.730 1 15473.413 14992.629
frequency modulation -3.158 -4.709 -1.594 1.002 2625.420 5199.726
amplitude modulation -2.724 -4.204 -1.227 1.001 3167.664 5593.286
frequency 3.164 2.902 3.431 1.002 2763.472 5433.162
duration -0.849 -2.367 0.687 1.001 2700.931 5700.373
frequency:habitat structure -0.639 -0.785 -0.493 1 25778.804 13335.766
habitat structure:duration -0.399 -1.241 0.427 1 20933.514 14843.771
habitat structure:amplitude modulation -1.576 -2.414 -0.739 1 21508.820 14345.505
habitat structure:frequency modulation -0.128 -0.967 0.711 1 22068.512 14167.521

3.3.1 Conditional plots

Code
bs <- 10
ts <- 4

ggf <- plot(conditional_effects(mod, "frequency:habitat.structure"),
    plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Open",
    "Closed")) + theme_classic(base_size = bs) + labs(x = "Frequency",
    y = "Degradation (PC1)", color = "Habitat") + scale_fill_discrete(guide = "none") +
    theme(legend.position = c(0.2, 0.95)) + scale_x_continuous(breaks = seq(-5,
    5, 2.5), labels = round(seq(-5, 5, 2.5) + mean(read.csv("./data/processed/tlalpan_degradation_metrics_v01.csv")$frequency),
    1)) + geom_text(x = 0, y = 70, label = "*", size = ts)

ylim <- c(20, 80)

ggfm <- plot(conditional_effects(mod, "habitat.structure:frequency.modulation"),
    plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Tonal",
    "FM")) + theme_classic(base_size = bs) + labs(x = "Habitat", y = "",
    color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
    0.5)) + ylim(ylim)

diff <- 0.1


# add contrasts
for (i in 1:2) ggfm <- ggfm + geom_signif(y_position = c(65), xmin = c(i -
    diff), xmax = c(i + diff), annotation = c("*"), tip_length = 0.02,
    textsize = ts, size = 0.4, color = "gray40")

ggfm <- ggfm + geom_signif(y_position = c(72), xmin = 1, xmax = 2,
    annotation = c("ns"), tip_length = 0.02, textsize = ts, size = 0.4,
    color = "gray40")

ggam <- plot(conditional_effects(mod, "habitat.structure:amplitude.modulation"),
    plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Flat",
    "AM")) + theme_classic(base_size = bs) + labs(x = "Habitat", y = "Degradation (PC1)",
    color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
    0.5)) + ylim(ylim)

# add contrasts
for (i in 1:2) ggam <- ggam + geom_signif(y_position = c(65), xmin = c(i -
    diff), xmax = c(i + diff), annotation = c("*"), tip_length = 0.02,
    textsize = ts, size = 0.4, color = "gray40")

ggam <- ggam + geom_signif(y_position = c(72), xmin = 1, xmax = 2,
    annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4,
    color = "gray40")


ggdur <- plot(conditional_effects(mod, "habitat.structure:duration"),
    plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Short",
    "Long")) + theme_classic(base_size = bs) + labs(x = "Habitat",
    y = "", color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
    0.5)) + ylim(ylim)

# add contrasts
for (i in 1:2) ggdur <- ggdur + geom_signif(y_position = c(65), xmin = c(i -
    diff), xmax = c(i + diff), annotation = c("ns"), tip_length = 0.02,
    textsize = ts, size = 0.4, color = "gray40")

ggdur <- ggdur + geom_signif(y_position = c(72), xmin = 1, xmax = 2,
    annotation = c("ns"), tip_length = 0.02, textsize = 4, size = 0.4,
    color = "gray40")


plot_grid(ggf, ggfm, ggam, ggdur)

3.4 Tail-to-signal ratio

Code
mod <- readRDS("./data/processed/regression_model_tail_to_signal_ratio.RDS")

extended_summary(fit = mod, n.posterior = 1000, fill = viridis(10)[7],
    trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE,
    gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam",
        "frequency.modulationfm", "durationlong"), gsub.replacement = c("",
        "habitat structure", "amplitude modulation", "frequency modulation",
        "duration"), effects = effects, print.name = FALSE)
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 6) Intercept-student_t(3, -6.3, 8.7) sd-cauchy(0, 4) sigma-student_t(3, 0, 8.7) simo-dirichlet(1) tail.to.signal.ratio ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates) 10000 4 1 5000 16 0 2380.891 4713.849 1732769608
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
habitat structure 3.253 2.897 3.605 1 13266.787 14830.753
frequency modulation -0.818 -1.436 -0.211 1.001 2503.012 5416.822
amplitude modulation 0.822 0.213 1.450 1.001 2720.556 5229.965
frequency 0.303 0.196 0.409 1.002 2749.731 4929.226
duration -0.426 -1.029 0.182 1.001 2380.891 4713.849
frequency:habitat structure 0.020 -0.041 0.082 1 25850.405 15259.070
habitat structure:duration -0.053 -0.412 0.301 1 18539.674 15164.141
habitat structure:amplitude modulation -0.385 -0.741 -0.030 1 17231.026 14672.205
habitat structure:frequency modulation 0.540 0.183 0.898 1 20156.979 15407.009

3.4.1 Conditional plots

Code
bs <- 10
ts <- 4

ggf <- plot(conditional_effects(mod, "frequency:habitat.structure"),
    plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Open",
    "Closed")) + theme_classic(base_size = bs) + labs(x = "Frequency",
    y = "Degradation (PC1)", color = "Habitat") + scale_fill_discrete(guide = "none") +
    theme(legend.position = c(0.2, 0.95)) + scale_x_continuous(breaks = seq(-5,
    5, 2.5), labels = round(seq(-5, 5, 2.5) + mean(read.csv("./data/processed/tlalpan_degradation_metrics_v01.csv")$frequency),
    1)) + geom_text(x = 0, y = 70, label = "*", size = ts)

ylim <- c(-27, -5)

ggfm <- plot(conditional_effects(mod, "habitat.structure:frequency.modulation"),
    plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Tonal",
    "FM")) + theme_classic(base_size = bs) + labs(x = "Habitat", y = "",
    color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
    0.5)) + ylim(ylim)

diff <- 0.1


# add contrasts
for (i in 1:2) ggfm <- ggfm + geom_signif(y_position = c(-12), xmin = c(i -
    diff), xmax = c(i + diff), annotation = c("*"), tip_length = 0.02,
    textsize = ts, size = 0.4, color = "gray40")

ggfm <- ggfm + geom_signif(y_position = c(-9), xmin = 1, xmax = 2,
    annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4,
    color = "gray40")

ggam <- plot(conditional_effects(mod, "habitat.structure:amplitude.modulation"),
    plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Flat",
    "AM")) + theme_classic(base_size = bs) + labs(x = "Habitat", y = "Degradation (PC1)",
    color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
    0.5)) + ylim(ylim)

# add contrasts
for (i in 1:2) ggam <- ggam + geom_signif(y_position = c(-12), xmin = c(i -
    diff), xmax = c(i + diff), annotation = c("*"), tip_length = 0.02,
    textsize = ts, size = 0.4, color = "gray40")

ggam <- ggam + geom_signif(y_position = c(-9), xmin = 1, xmax = 2,
    annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4,
    color = "gray40")


ggdur <- plot(conditional_effects(mod, "habitat.structure:duration"),
    plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Short",
    "Long")) + theme_classic(base_size = bs) + labs(x = "Habitat",
    y = "", color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
    0.5)) + ylim(ylim)

# add contrasts
for (i in 1:2) ggdur <- ggdur + geom_signif(y_position = c(-12), xmin = c(i -
    diff), xmax = c(i + diff), annotation = c("ns"), tip_length = 0.02,
    textsize = ts, size = 0.4, color = "gray40")

ggdur <- ggdur + geom_signif(y_position = c(-9), xmin = 1, xmax = 2,
    annotation = c("ns"), tip_length = 0.02, textsize = 4, size = 0.4,
    color = "gray40")


plot_grid(ggf, ggfm, ggam, ggdur)

3.5 Combined results

Code
pc1 <- extended_summary(read.file = "./data/processed/regression_model_pc1.RDS",
    n.posterior = 1000, fill = viridis(10)[7], trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects, print.name = FALSE, trace = FALSE, return = TRUE)

gg_pc1 <- pc1$plot + labs(x = "Degradation (PC1)")

br <- extended_summary(read.file = "./data/processed/regression_model_blur_ratio.RDS",
    n.posterior = 1000, fill = viridis(10)[7], trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects, print.name = FALSE, trace = FALSE, return = TRUE)

gg_br <- br$plot + labs(x = "Blur ratio", y = "") + theme(axis.text.y = element_blank())

ea <- extended_summary(read.file = "./data/processed/regression_model_excess_attenuation.RDS",
    n.posterior = 1000, fill = viridis(10)[7], trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects, print.name = FALSE, trace = FALSE, return = TRUE)

gg_ea <- ea$plot + labs(x = "Excess attenuation", y = "") + theme(axis.text.y = element_blank())

tsr <- extended_summary(read.file = "./data/processed/regression_model_tail_to_signal_ratio.RDS",
    n.posterior = 1000, fill = viridis(10)[7], trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects, print.name = FALSE, trace = FALSE, return = TRUE)

gg_tsr <- tsr$plot + labs(x = "Tail-to-signal ratio", y = "") + theme(axis.text.y = element_blank())

plot_grid(gg_pc1, gg_br, gg_ea, gg_tsr, nrow = 1, rel_widths = c(6,
    2, 2, 2))

Code
estimates <- data.frame(pc1 = pc1$coef_table$Estimate, br = br$coef_table$Estimate,
    ea = ea$coef_table$Estimate, tsr = tsr$coef_table$Estimate)

names(estimates) <- c("Degradation (PC1)", "Blur ratio", "Excess attenuation",
    "Tail-to-signal ratio")

# make estimates relative to maximum estimate in data
rel_estimates <- data.frame(lapply(estimates, function(x) x/max(abs(x)) *
    0.8))

names(rel_estimates) <- c("Degradation (PC1)", "Blur ratio", "Excess attenuation",
    "Tail-to-signal ratio")

# corrplot(as.matrix(estimates), method = 'ellipse', )

signif <- data.frame(pc1 = pc1$coef_table$`l-95% CI` * pc1$coef_table$`u-95% CI` >
    0, br = br$coef_table$`l-95% CI` * br$coef_table$`u-95% CI` >
    0, ea = ea$coef_table$`l-95% CI` * ea$coef_table$`u-95% CI` >
    0, tsr = tsr$coef_table$`l-95% CI` * tsr$coef_table$`u-95% CI` >
    0)

names(signif) <- c("Degradation (PC1)", "Blur ratio", "Excess attenuation",
    "Tail-to-signal ratio")

estimates <- cbind(rownames(pc1$coef_table), stack(estimates)[, 1],
    stack(rel_estimates)[, 1], stack(signif))

names(estimates) <- c("predictor", "est", "relavite.est", "sig", "response")
estimates$relavite.est <- ifelse(estimates$sig, estimates$relavite.est,
    0)

estimates$response <- factor(estimates$response, levels = rev(c("Degradation (PC1)",
    "Blur ratio", "Excess attenuation", "Tail-to-signal ratio")))

estimates$predictor <- factor(estimates$predictor, levels = c("habitat structure",
    "frequency modulation", "amplitude modulation", "frequency", "duration",
    "frequency:habitat structure", "habitat structure:duration", "habitat structure:amplitude modulation",
    "habitat structure:frequency modulation"))

ggplot(estimates, aes(x = predictor, y = response, fill = relavite.est)) +
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + geom_tile()
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + +
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + coord_equal()
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + +
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + #
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + use
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + blue
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + and
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + red
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + palette
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + scale_fill_gradient2(low
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + =
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + '#F36c68',
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + high
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + =
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + '#175493',
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + guide
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + =
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + 'none')
    geom_tile() + coord_equal() + # use blue and red palette  scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + +
scale_fill_gradient2(low = viridis(10)[3], high = viridis(10)[7],
    guide = "none") + geom_text(aes(label = round(est, 2)), color = "black") +
    labs(x = "", y = "") + theme_classic() + theme(axis.text.x = element_text(color = "black",
    size = 11, angle = 30, vjust = 0.8, hjust = 0.8))

 

4 Takeaways

  • Habitat structure seems to drive most of the observed degradation
  • Frequency modulation and amplitude modulation were the signal structural features that more strongly affected transmission
  • Signal features interact with habitat structure in diverse ways, sometimes increasing and sometimes decreasing degradation
  • Degradation metrics are robust to variation in signal duration

 

5 Posterior predictive checks

Code
model_list <- c(pc1 = "./data/processed/regression_model_pc1.RDS",
    blur_ratio = "./data/processed/regression_model_blur_ratio.RDS",
    excess_attenuation = "./data/processed/regression_model_excess_attenuation.RDS",
    tail_to_signal_ratio = "./data/processed/regression_model_tail_to_signal_ratio.RDS")

ndraws <- 20

for (i in seq_len(length(model_list))) {
    cat("\n\n## ", names(model_list)[i], "\n\n")

    fit <- readRDS(model_list[[i]])
    ppc_dens <- pp_check(fit, ndraws = ndraws, type = "dens_overlay_grouped",
        group = "habitat.structure")  # shows dens_overlay plot by default
    pp_mean <- pp_check(fit, type = "stat_grouped", stat = "mean",
        group = "habitat.structure", ndraws = ndraws)
    pp_scat <- pp_check(fit, type = "scatter_avg_grouped", group = "habitat.structure",
        ndraws = ndraws)
    pp_stat2 <- pp_check(fit, type = "stat_2d", ndraws = ndraws)
    pp_loo <- pp_check(fit, type = "loo_pit_qq", ndraws = ndraws)
    pp_error <- pp_check(fit, type = "error_scatter_avg_grouped",
        ndraws = ndraws, group = "habitat.structure")
    plot_l <- list(ppc_dens, pp_mean, pp_scat, pp_error, pp_stat2,
        pp_loo)

    plot_l <- lapply(plot_l, function(x) x + scale_color_viridis_d(begin = 0.1,
        end = 0.8, alpha = 0.5) + scale_fill_viridis_d(begin = 0.1,
        end = 0.8, alpha = 0.5) + theme_classic())

    print(plot_grid(plotlist = plot_l, ncol = 2))
}

5.1 pc1

5.2 blur_ratio

5.3 excess_attenuation

5.4 tail_to_signal_ratio


 

Session information

R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3

locale:
 [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=pt_BR.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggsignif_0.6.2     emmeans_1.8.1-1    brmsish_1.0.0      cmdstanr_0.4.0    
 [5] cowplot_1.1.1      ggdist_3.2.0       brms_2.18.0        Rcpp_1.0.10       
 [9] corrplot_0.90      viridis_0.6.2      viridisLite_0.4.1  tidyr_1.1.3       
[13] ggplot2_3.4.0      warbleR_1.1.28     NatureSounds_1.0.4 seewave_2.2.0     
[17] tuneR_1.4.4        rprojroot_2.0.3    formatR_1.11       knitr_1.43        
[21] kableExtra_1.3.4   remotes_2.4.2     

loaded via a namespace (and not attached):
  [1] backports_1.4.1      systemfonts_1.0.4    plyr_1.8.7          
  [4] igraph_1.4.3         splines_4.1.0        crosstalk_1.2.0     
  [7] TH.data_1.1-0        rstantools_2.2.0     inline_0.3.19       
 [10] digest_0.6.31        htmltools_0.5.5      fansi_1.0.4         
 [13] magrittr_2.0.3       checkmate_2.1.0      RcppParallel_5.1.5  
 [16] matrixStats_0.62.0   xts_0.12.2           sandwich_3.0-1      
 [19] svglite_2.1.0        prettyunits_1.1.1    colorspace_2.0-3    
 [22] signal_0.7-7         rvest_1.0.3          xfun_0.39           
 [25] dplyr_1.0.10         callr_3.7.3          crayon_1.5.2        
 [28] RCurl_1.98-1.12      jsonlite_1.8.5       lme4_1.1-27.1       
 [31] ape_5.6-2            survival_3.2-11      zoo_1.8-11          
 [34] glue_1.6.2           gtable_0.3.1         webshot_0.5.4       
 [37] distributional_0.3.1 pkgbuild_1.4.0       rstan_2.21.7        
 [40] abind_1.4-5          scales_1.2.1         mvtnorm_1.1-3       
 [43] DBI_1.1.3            miniUI_0.1.1.1       dtw_1.23-1          
 [46] xtable_1.8-4         proxy_0.4-27         stats4_4.1.0        
 [49] StanHeaders_2.21.0-7 DT_0.26              htmlwidgets_1.5.4   
 [52] httr_1.4.4           threejs_0.3.3        posterior_1.3.1     
 [55] ellipsis_0.3.2       pkgconfig_2.0.3      loo_2.4.1.9000      
 [58] farver_2.1.1         utf8_1.2.3           labeling_0.4.2      
 [61] tidyselect_1.2.0     rlang_1.1.1          reshape2_1.4.4      
 [64] later_1.3.0          munsell_0.5.0        tools_4.1.0         
 [67] cli_3.6.1            generics_0.1.3       ggridges_0.5.4      
 [70] evaluate_0.21        stringr_1.5.0        fastmap_1.1.1       
 [73] yaml_2.3.7           processx_3.8.1       purrr_1.0.0         
 [76] pbapply_1.7-0        nlme_3.1-152         mime_0.12           
 [79] projpred_2.0.2       xml2_1.3.3           brio_1.1.3          
 [82] compiler_4.1.0       bayesplot_1.9.0      shinythemes_1.2.0   
 [85] rstudioapi_0.14      gamm4_0.2-6          testthat_3.1.8      
 [88] tibble_3.2.1         stringi_1.7.12       ps_1.7.5            
 [91] Brobdingnag_1.2-9    lattice_0.20-44      Matrix_1.5-1        
 [94] nloptr_1.2.2.2       markdown_1.3         shinyjs_2.1.0       
 [97] fftw_1.0-7           tensorA_0.36.2       vctrs_0.6.2         
[100] pillar_1.9.0         lifecycle_1.0.3      bridgesampling_1.1-2
[103] estimability_1.4.1   bitops_1.0-7         httpuv_1.6.6        
[106] R6_2.5.1             promises_1.2.0.1     gridExtra_2.3       
[109] codetools_0.2-18     boot_1.3-28          colourpicker_1.2.0  
[112] MASS_7.3-54          gtools_3.9.3         assertthat_0.2.1    
[115] rjson_0.2.21         withr_2.5.0          shinystan_2.6.0     
[118] multcomp_1.4-17      mgcv_1.8-36          parallel_4.1.0      
[121] grid_4.1.0           minqa_1.2.4          coda_0.19-4         
[124] rmarkdown_2.20       shiny_1.7.3          base64enc_0.1-3     
[127] dygraphs_1.1.1.6